In [1]:
%matplotlib inline
%config InlineBackend.figure_format='retina'
# %config InlineBackend.figure_format='svg'
import os
import time
from glob import glob
import numpy as np
brick = 23
STARFISH = os.getenv("STARFISH")
isoc_dir = "b23ir_isoc"
lib_dir = "b23ir_lib"
synth_dir = "b23ir_synth"
fit_dir = "b23ir_fit"
wfc3_bands = ['F110W', 'F160W']
In [2]:
from astropy.coordinates import Distance
import astropy.units as u
from padova import AgeGridRequest, IsochroneRequest
from starfisher import LibraryBuilder
z_grid = [0.012, 0.015, 0.019, 0.024, 0.03]
delta_gyr = 0.5
late_ages = np.log10(np.arange(1e9 + delta_gyr, 13e9, delta_gyr * 1e9))
print late_ages
if not os.path.exists(os.path.join(STARFISH, isoc_dir)):
for z in z_grid:
# Young ages in log-grid
r = AgeGridRequest(z,
min_log_age=6.6,
max_log_age=9.,
delta_log_age=0.5,
phot='wfc3', photsys_version='odfnew')
for isoc in r.isochrone_set:
isoc.export_for_starfish(os.path.join(STARFISH, isoc_dir),
bands=wfc3_bands)
# Old ages in linear grid
for logage in late_ages:
r = IsochroneRequest(z, logage,
phot='wfc3', photsys_version='odfnew')
r.isochrone.export_for_starfish(os.path.join(STARFISH, isoc_dir),
bands=wfc3_bands)
d = Distance(785 * u.kpc)
builder = LibraryBuilder(isoc_dir, lib_dir,
nmag=len(wfc3_bands),
dmod=d.distmod.value,
iverb=3)
if not os.path.exists(builder.full_isofile_path):
builder.install()
In [3]:
from collections import namedtuple
from starfisher import Lockfile
from starfisher import Synth
from starfisher import ExtinctionDistribution
from starfisher import ExtantCrowdingTable
from starfisher import ColorPlane
from m31hst.phatast import PhatAstTable
if not os.path.exists(os.path.join(STARFISH, synth_dir)):
os.makedirs(os.path.join(STARFISH, synth_dir))
# No binning in our lockfile
lockfile = Lockfile(builder.read_isofile(), synth_dir, unbinned=True)
# No extinction, yet
young_av = ExtinctionDistribution()
old_av = ExtinctionDistribution()
rel_extinction = np.ones(len(wfc3_bands), dtype=float)
for av in (young_av, old_av):
av.set_uniform(0.)
# Use PHAT AST from the outer field
crowd_path = os.path.join(synth_dir, "crowding.dat")
full_crowd_path = os.path.join(STARFISH, crowd_path)
tbl = PhatAstTable()
tbl.write_crowdfile_for_field(full_crowd_path, 0,
bands=('f110w', 'f160w'))
crowd = ExtantCrowdingTable(crowd_path)
# Define CMD planes
Lim = namedtuple('Lim', 'x y')
ir_lim = Lim(x=(-1., 2.), y=(25., 16.))
ir_cmd = ColorPlane((wfc3_bands.index('F110W'),
wfc3_bands.index('F160W')),
wfc3_bands.index('F160W'),
ir_lim.x,
(min(ir_lim.y), max(ir_lim.y)),
28.,
suffix='f110f160',
x_label=r'$\mathrm{F110W}-\mathrm{F160W}$',
y_label=r'$\mathrm{F110W}$',
dpix=0.1)
colour_planes = [ir_cmd]
synth = Synth(synth_dir, builder, lockfile, crowd,
rel_extinction,
young_extinction=young_av,
old_extinction=old_av,
planes=colour_planes,
mass_span=(0.08, 150.),
nstars=10000000)
if len(glob(os.path.join(STARFISH, synth_dir, "z*"))) == 0:
synth.run_synth(n_cpu=4)
synth.plot_all_hess(os.path.join(STARFISH, synth_dir, 'hess'))
In [4]:
from astropy.table import Table
from m31hst import phat_v2_phot_path
if not os.path.exists(os.path.join(STARFISH, fit_dir)):
os.makedirs(os.path.join(STARFISH, fit_dir))
data_root = os.path.join(fit_dir, "b23ir.")
full_data_path = os.path.join(STARFISH, '{0}f110f160'.format(data_root))
brick_table = Table.read(phat_v2_phot_path(brick), format='fits')
# Only use stars within the fitting box
c = brick_table['f110w_vega'] - brick_table['f160w_vega']
m = brick_table['f160w_vega']
sel = np.where((c > min(ir_lim.x)) & (c < max(ir_lim.x)) &
(m > min(ir_lim.y)) & (m < max(ir_lim.y)))[0]
brick_table = brick_table[sel]
print("Fitting {0:d} stars".format(len(brick_table)))
if not os.path.exists(full_data_path):
phot_dtype = np.dtype([('x', np.float), ('y', np.float)])
photdata = np.empty(len(brick_table), dtype=phot_dtype)
photdata['x'][:] = brick_table['f110w_vega'] - brick_table['f160w_vega']
photdata['y'][:] = brick_table['f160w_vega']
np.savetxt(full_data_path, photdata, delimiter=' ', fmt='%.4f')
In [5]:
import matplotlib as mpl
import matplotlib.pyplot as plt
from androcmd.plot import contour_hess
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
contour_hess(ax, brick_table['f110w_vega'] - brick_table['f160w_vega'],
brick_table['f160w_vega'], ir_lim.x, (max(ir_lim.y), min(ir_lim.y)),
plot_args={'ms': 3})
ir_cmd.plot_mask(ax)
ax.set_xlabel(r'$\mathrm{F110W}-\mathrm{F160W}$')
ax.set_ylabel(r'$\mathrm{F110W}$')
ax.set_xlim(ir_lim.x)
ax.set_ylim(ir_lim.y)
fig.show()
In [6]:
from starfisher import SFH, Mask
mask = Mask(colour_planes)
sfh = SFH(data_root, synth, mask, fit_dir)
if not os.path.exists(sfh.full_outfile_path):
sfh.run_sfh()
sfh_table = sfh.solution_table()
print(sfh_table)
In [9]:
from starfisher.sfhplot import SFHCirclePlot
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
cp = SFHCirclePlot(sfh_table)
cp.plot_in_ax(ax, max_area=800)
for logage in np.log10(np.arange(1, 13, 1) * 1e9):
ax.axvline(logage, c='0.8', zorder=-1)
ax.set_ylim(-0.5, 0.4)
fig.show()
In [8]:
import cubehelix
cmapper = lambda: cubehelix.cmap(startHue=240,endHue=-300,minSat=1,maxSat=2.5,minLight=.3,maxLight=.8,gamma=.9)
from starfisher.sfhplot import ChiTriptykPlot
fig = plt.figure(figsize=(7, 5))
ctp = ChiTriptykPlot(sfh.full_chi_path, 1, ir_cmd.x_span, ir_cmd.y_span,
ir_cmd.dpix, ir_cmd.x_label, ir_cmd.y_label,
flipy=True)
ax_obs, ax_mod, ax_chi = ctp.setup_axes(fig)
ctp.plot_obs_in_ax(ax_obs, cmap=cmapper())
ctp.plot_mod_in_ax(ax_mod, cmap=cmapper())
ctp.plot_chi_in_ax(ax_chi, cmap=cubehelix.cmap())
ax_obs.text(0.0, 1.01, "Observed", transform=ax_obs.transAxes, size=8, ha='left')
ax_mod.text(0.0, 1.01, "Model", transform=ax_mod.transAxes, size=8, ha='left')
ax_chi.text(0.0, 1.01, r"$\log \chi^2$", transform=ax_chi.transAxes, size=8, ha='left')
fig.show()
In [ ]: